All samples

design <- model.matrix(~ 0 + Group, data = salmon@phenoData)
colnames(design) <- gsub(pattern = "Group", replacement = "", x = colnames(design))

y <- DGEList(counts = salmon@gene.counts)
keep <- filterByExpr(y, design, min.count = 20)
y <- y[keep, ]

dds <- calcNormFactors(y)
v <- voom(dds, design = design)

contrast.matrix <- makeContrasts(MSUS - CTRL, levels = design)
fit <- lmFit(v)

fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)

dea <- as.data.frame(topTable(fit2, coef = 1, number = Inf))
dea <- data.frame(Genes = rownames(dea), dea, stringsAsFactors = F)
dea_sig <- dea[abs(dea$logFC) >= 1 & dea$adj.P.Val <= 0.05, ]

Removing L001

design_l1 <- model.matrix(~ 0 + Group, data = salmon@phenoData[-1,])
colnames(design_l1) <- gsub(pattern = "Group", replacement = "", x = colnames(design_l1))

y_l1 <- DGEList(counts = salmon@gene.counts[,-1])
keep_l1 <- filterByExpr(y_l1, design_l1, min.count = 20)
y_l1 <- y_l1[keep_l1, ]

dds_l1 <- calcNormFactors(y_l1)
v_l1 <- voom(dds_l1, design = design_l1)

contrast.matrix_l1 <- makeContrasts(MSUS - CTRL, levels = design_l1)
fit_l1 <- lmFit(v_l1)

fit2_l1 <- contrasts.fit(fit_l1, contrast.matrix_l1)
fit2_l1 <- eBayes(fit2_l1)

dea_l1 <- as.data.frame(topTable(fit2_l1, coef = 1, number = Inf))
dea_l1 <- data.frame(Genes = rownames(dea_l1), dea_l1, stringsAsFactors = F)
dea_sig_l1 <- dea_l1[abs(dea_l1$logFC) >= 1 & dea_l1$adj.P.Val <= 0.05, ]

Save RData files

save(dea.list,
  file = "./output/dea_mPFC_lit.DEA.RData", compress = T,
  compression_level = 3
)

save(normCounts,
  file = "./output/normCounts_voom_mPFC_lit.RData", compress = T,
  compression_level = 3
)

save(voomEList,
  file = "./output/voom_EList_mPFC_lit.RData", compress = T,
  compression_level = 3
)

save(dea,
  file = "./output/limma_mPFC_lit.RData", compress = T,
  compression_level = 3
)

writexl::write_xlsx(x = dea, path = "output/dea_results.xlsx", col_names = T, format_headers = T)

All samples

res <- dea
res$threshold <- as.factor(res$adj.P.Val < 0.1)

p1 <- ggplot(data = res, aes(
  x = res$AveExpr,
  y = res$logFC,
  colour = threshold
)) +
  geom_point(alpha = 0.5, size = 1.8) +
  geom_hline(aes(yintercept = 0), colour = "blue", size = 1) +
  ylim(c(
    -ceiling(max(abs(res$logFC))),
    ceiling(max(abs(res$logFC)))
  )) +
  ggtitle("MSUS vs CTRL") +
  labs(subtitle = "loess fit") +
  xlab("Mean expression") +
  ylab("Log2 Fold Change") +
  theme_bw() +
  theme(
    plot.title = element_text(face = "bold", size = 20, hjust = 0.5),
    plot.subtitle = element_text(size = 16, hjust = 0.5, face = "italic", color = "blue"),
    axis.title = element_text(face = "bold", size = 15),
    axis.text = element_text(face = "bold", size = 12),
    legend.title = element_text(face = "bold", size = 15),
    legend.text = element_text(size = 14),
    panel.border = element_rect(colour = "black", fill=NA, size=2)
  ) +
  scale_color_manual(name = "p.adjusted < 0.1", values=brewer.pal(n = 3, name = "Set1")[1:2])+ NULL
  stat_smooth(se = FALSE, method = "loess", color = "red", formula = y ~ x, size = 1)

## geom_smooth: se = FALSE, na.rm = FALSE, orientation = NA
## stat_smooth: method = loess, formula = y ~ x, se = FALSE, n = 80, fullrange = FALSE, level = 0.95, na.rm = FALSE, orientation = NA, method.args = list(), span = 0.75
## position_identity

Removing L0001

res_l1 <- dea_l1
res_l1$threshold <- as.factor(res_l1$adj.P.Val < 0.1)

p2 <- ggplot(data = res_l1, aes(
  x = res_l1$AveExpr,
  y = res_l1$logFC,
  colour = threshold
)) +
  geom_point(alpha = 0.8, size = 3) +
  geom_hline(aes(yintercept = 0), colour = "blue", size = 1) +
  ylim(c(
    -ceiling(max(abs(res_l1$logFC))),
    ceiling(max(abs(res_l1$logFC)))
  )) +
  ggtitle("MSUS vs CTRL (no L1)") +
  labs(subtitle = "loess fit") +
  xlab("Mean expression") +
  ylab("Log2 Fold Change") +
  theme_bw() +
  theme(
    plot.title = element_text(face = "bold", size = 20, hjust = 0.5),
    plot.subtitle = element_text(size = 16, hjust = 0.5, face = "italic", color = "blue"),
    axis.title = element_text(face = "bold", size = 15),
    axis.text = element_text(face = "bold", size = 12),
    legend.title = element_text(face = "bold", size = 15),
    legend.text = element_text(size = 14),
    panel.border = element_rect(colour = "black", fill=NA, size=2)
  ) +
  scale_color_manual(name = "p.adjusted < 0.1", values=brewer.pal(n = 3, name = "Set1")[1:2])+ NULL
  stat_smooth(se = FALSE, method = "loess", color = "red", formula = y ~ x, size = 1)

## geom_smooth: se = FALSE, na.rm = FALSE, orientation = NA
## stat_smooth: method = loess, formula = y ~ x, se = FALSE, n = 80, fullrange = FALSE, level = 0.95, na.rm = FALSE, orientation = NA, method.args = list(), span = 0.75
## position_identity

p1 | p2

## Warning: Use of `res$AveExpr` is discouraged. Use `AveExpr` instead.

## Warning: Use of `res$logFC` is discouraged. Use `logFC` instead.

## Warning: Use of `res_l1$AveExpr` is discouraged. Use `AveExpr` instead.

## Warning: Use of `res_l1$logFC` is discouraged. Use `logFC` instead.

All samples

xvar <- apply((data$cpm + 1), 1, var)
genes500 <- head(sort(xvar, decreasing = TRUE), n = 500)
x500 <- data$cpm[rownames(data$cpm) %in% names(genes500), ]

pheatmap(x500, scale = "row", show_rownames = F, col = viridis(100), show_colnames = F, annotation_col = data$pData[,2:3], main = "Clustering of samples based on top 500 variable genes")

Removing L001

xvar <- apply((data$cpm[,-1] + 1), 1, var)
genes500 <- head(sort(xvar, decreasing = TRUE), n = 500)
x500 <- data$cpm[rownames(data$cpm) %in% names(genes500), ]

pheatmap(x500, scale = "row", show_rownames = F, col = viridis(100), show_colnames = F, annotation_col = data$pData[,2:3], main = "Clustering of samples based on top 500 variable genes (no L1)")

All samples

pca <- prcomp(t(data$cpm))
# ggplotly(
  pc1 <- autoplot(pca,
  data = data$pData, colour = "Group",
  frame = TRUE, frame.type = "norm", size = 5
) +
  ggtitle("") + 
    theme_bw() +
    scale_color_manual(values=brewer.pal(n = 3, name = "Set1")[1:2]) +
    theme(axis.text.x = element_text(size = 20, colour = "black"),
        axis.text.y = element_text(size = 20, colour = "black"),  
        axis.title.x = element_text(size = 20, colour = "black"),
        axis.title.y = element_text(size = 20, colour = "black"),
        plot.title = element_text(size = 25, face = "bold", colour = "black"),
        legend.text=element_text(size=20, colour = "black"),
        legend.title=element_text(size=20, face = "bold"),
        panel.border = element_rect(colour = "black", fill=NA, size=2))
    
  
  pc2 <- autoplot(pca,
  data = data$pData, colour = "Lane",
  frame = TRUE, frame.type = "norm", size = 5
) +
  ggtitle("") + 
    theme_bw()+
    scale_color_manual(values=brewer.pal(n = 3, name = "Set1")[1:2]) +
    theme(axis.text.x = element_text(size = 20),
        axis.text.y = element_text(size = 20),  
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20),
        plot.title = element_text(size = 25, face = "bold"),
        legend.text=element_text(size=20),
        legend.title=element_text(size=20, face = "bold"),
        panel.border = element_rect(colour = "black", fill=NA, size=2))
  
  
  pca1 <- prcomp(t(data$cpm[,-1]))
# ggplotly(
  pc3 <- autoplot(pca1,
  data = data$pData[-1,], colour = "Group",
  frame = TRUE, frame.type = "norm", size = 5
) +
  ggtitle("") + 
    theme_bw() +
    scale_color_manual(values=brewer.pal(n = 3, name = "Set1")[1:2]) +
    theme(axis.text.x = element_text(size = 20, colour = "black"),
        axis.text.y = element_text(size = 20, colour = "black"),  
        axis.title.x = element_text(size = 20, colour = "black"),
        axis.title.y = element_text(size = 20, colour = "black"),
        plot.title = element_text(size = 25, face = "bold", colour = "black"),
        legend.text=element_text(size=20, colour = "black"),
        legend.title=element_text(size=20, face = "bold"),
        panel.border = element_rect(colour = "black", fill=NA, size=2))

pc1|pc2

## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse

pc3

## Too few points to calculate an ellipse
## Too few points to calculate an ellipse

All samples

pca1 <- plgINS::plPCA(x = data$cpm, samples_data = data$pData, colorBy = "Group", points.size = 20)

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:IRanges':
## 
##     slice

## The following object is masked from 'package:S4Vectors':
## 
##     rename

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

## Warning: `arrange_()` is deprecated as of dplyr 0.7.0.
## Please use `arrange()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.

## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels

## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels

## Warning: Didn't find a colorbar to modify.

Remove L1

pca2 <- plgINS::plPCA(x = data$cpm[,-1], samples_data = data$pData[-1,], colorBy = "Group", points.size = 20)

## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels

## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels

## Warning: Didn't find a colorbar to modify.

subplot(pca1, pca2)